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ABSTRACT 

The Crab Nebula was detected with the Milagro experiment at a statistical signif- 
icance of 17 standard deviations over the lifetime of the experiment. The experiment 
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was sensitive to approximately 100 GeV -100 TeV gamma ray air showers by observing 
the particle footprint reaching the ground. The fraction of detectors recording signals 
from photons at the ground is a suitable proxy for the energy of the primary particle 
and has been used to measure the photon energy spectrum of the Crab Nebula between 
and ~100 TeV. The TeV emission is believed to be caused by inverse- Compton 
up-scattering scattering of ambient photons by an energetic electron population. The 
location of a TeV steepening or cutoff in the energy spectrum reveals important de- 
tails about the underlying electron population. We describe the experiment and the 
technique for distinguishing gamma-ray events from the much more-abundant hadronic 
events. We describe the calculation of the significance of the excess from the Crab and 
how the energy spectrum is fit. The excess from the Crab is fit to a function of the form 

^{lo.a.Ecut) = lo(^) 

where the flux /q, spectral index a and exponential cutoff energy £^cut ^3:e allowed to 
vary and Eq is chosen to de-correlate the fit parameters. The energy spectrum, including 
the statistical errors from the fit, obtained using the simple power law hypothesis, that 
is Ecut “ oc, for data between September 2005 and March 2008 is: 

% = (6.5 ± 0.4) X 10-'^^(E/10 TeV)-3 i±o.i^Pj^2 TeV)”^ 

between TeV and '^100 TeV. When a finite Ecut Is fit the result is 

^ = (2.5+;^;^) X 10 - 12 (£;/3TeV)-2-5±o.4gj^p(_£:/32+39 TeV)(cm2 secTeV)“^ 

The results are subject to an ~30% systematic uncertainty in the overall flux and 
an ^0.1 in the power law indie ies quoted. Uncertainty in the overall energy scale has 
been absorbed into these errors. 

Fixing the spectral index to values that have been measured below 1 TeV by I ACT 
experiments (2.4 to 2.6), the fit to the Milagro data suggests that Crab exhibits a 
spectral steepening or cutoff between about 20 to 40 TeV. 

Subject headings: gamma rays: observations — pulsars: general — pulsar wind nebulae 
— Crab Nebula 


1. Introduction 

The Crab supernova remnant is a luminous neai'by VHE gamma-ray source created by a 
1054 CE supernova observed on Earth by Chinese, Arab and native American astronomers. The 
optically luminous shell is easily visible from Earth. The Crab Nebula lies 2 kpc from the Earth and 
is powered by a 33 ms pulsar that injects relativistic electrons into the nebula. The central pulsar 
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and its surrounding nebula are among the most widely studied astronomical objects across the 
entire electromagnetic spectrum. A population of energetic electrons is created by the conversion 
of rotational kinetic energy of the neutron star and acceleration in the shock formed where this 
flow reaches the surrounding medium. These electrons can interact by the inverse-Compton process 
with the associated synchrotron photons to create the multi- TeV gamma rays that have been seen 
(Gaensler & Slane 2006). 

Despite the recent flares in 100 MeV - 100 GeV emission observed in AGILE (Tavani et al. 
2011) and the Fermi-LAT (Abdo et al. 2011) and in the TeV by ARGO-YBJ (Aielli et al. 2010), 
the TeV emission from the Crab is believed to be steady when observed over several months and is 
a standard reference source for comparison to other TeV instruments. Such comparisons are useful 
as a cross-calibration of the various ground telescopes. The Crab was first detected at TeV energies 
by the Whipple telescope in 1989 (Weekes et al. 1989). Imaging Atmospheric Cherenkov Telescopes 
(lACTs) (Aharonian et al. 2004, 2006; Celik 2008: Aliu et al. 2008) and extensive air-shower (EAS) 
ground arrays (Atkins et al. 2003; Amenomori et al. 2009; Bartoli et al. 2011) have been used to 
identify and measure the flux of the TeV emission from the Crab. Since the size of the Crab nebula 
is small compared to the point spread function of TeV gamma-ray detectors, the emission region 
appears point-like, and study of the Crab can serve as a calibration of an instrument’s pointing 
and angular resolution. 

The Milagro experiment (Atkins et al. 2003; Abdo et al. 2009) was a large water-Cherenkov 
detector sensitive to energetic secondary particles in the particle shower resulting when a high- 
energy gamma ray or cosmic ray strikes the atmosphere. The experiment was able to distinguish 
gamma-ray-induced showers from hadron-induced showers by measuring the penetrating component 
characteristic of hadronic particle showers. The experiment was sensitive to EAS resulting from 
primary gamma rays between 100 GeV and 100 TeV and has dynamic range to resolve gamma-ray 
energy spectra between about 1 and 100 TeV. The experiment operated nearly 24 hours a day and 
viewed the entire overhead sky. The detector was located at 106.68^W longitude, 35.88^N latitude 
in northern New Mexico at an altitude of 2630 m above sea level and operated from 2000 to 2008. 
The sensitive area of the detector comprised two parts; a central water reservoir and an “outrigger” 
array of water tanks, both instrumented with 8- inch hemispherical photomultiplier tubes (PMTs) 
manufactured by Hamamatsu (Model R5912). The central reservoir was operated alone from 2000 
to 2004 when the outrigger array was added. The central Milagro reservoir consisted of two PMT 
layers (Atkins et al. 2000) deployed in a 60x80x8-meter covered water reservoir. The outrigger 
array consisted of 175 water tanks each with a single PMT mounted at the top of the tank and 
observing downward into the water of the Twek- lined tank. The outrigger tanks were spread in an 
irregular pattern over an area of 200x200 meters around the central reservoir. The PMTs detected 
Cherenkov radiation produced in the water by the high-energy particles in an extensive air shower 
(EAS) that reached to ground level. The Milagro experiment has previously reported TeV emission 
from the Crab (Atkins et al. 2003) prior to the addition of the Milagro outriggers. With the greater 
sensitivity from the outriggers and additional exposure, measurement of an energy spectrum is 



possible. The direction of the primary particle in an EAS was estimated using the arrival time of 
the PMT signals. The angular resolution of Milagro, defined as the standard deviation, cr, of a 
two-dimensional Gaussian fit to the angular error distribution, varied between about 1,2 and 0.35 
degrees for the data presented here. The angular resolution is a function of both the size of the 
event, and the operational period of the detector. 

Since higher-energy primary particles result in characteristically larger events on the ground, we 
use a measure of size of the events on the ground to measure the spectrum of the Crab, constraining 
emission out to 100 TeV. In section 2, we describe the background estimation and the construction 
of "skymaps’b Section 3 describes the event energy estimation and spectral fitting. In section 4, 
we verify our energy reconstruction using cosmic-ray hadrons and finally compute the flux and the 
spectrum of the Crab from 1 to 100 TeV. 


2. Skymaps and Background Estimation 

From the reconstructed data, a skymap - a histogram of the sky containing the number 
of events originating from each location and associated errors - is formed. These skymaps are 
binned in units of 0.1 deg and cover the viewable sky. All events are recorded in the J2000 epoch. 
Each recorded skymap contains a signal map and a background map which contain the measured 
counts on the sky and the background expectation respectively. The skymaps are constructed in 
independent bins of energy parameter T which is defined below in Equation 3. 

The hadronic background flux is stable in time at TeV energies because TeV cosmic-rays 
originate from distant sources and propagate diffusively in Galactic magnetic fields. Therefore, 
the TeV hadronic background is not strongly affected by local variations such as solar activity. 
Instead, the rate and angular distribution of events is dominated by variations in the atmosphere 
and the detector. The background computation technique described below is intended to measure 
and correct for these changes. 

The panels of Figure 1 demonstrate an example of the background computation in a single 
declination band. We represent the background rate F(r, /i, S) as a function of sidereal time r 
and declination 5 and local hour angle h. To great precision, F(r, /i, 6) can be separated into two 
independent terms, R(r)-e{h, S), where R{r) is the all-sky event rate and e{h. S) is the local angular 
distribution of events. 

Even large changes in due to trigger threshold changes for example, lead to only small 

changes in the angular distribution of events on the sky, e(h,5). We exploit this feature of at- 
mospheric showers to compute the background B{a. d) in celestial coordinates right-ascension, a, 
and declination, 8. The technique begins with the definition of an integration duration. For most 
data, the integration duration is two hours but when looking at rare events (very hard cuts) the 
integration duration is 24 hours. acquire data for the integration duration and form the effi- 
ciency map eih>, 6). This efficiency map is a normalized probability density function which indicates 
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from where, in local detector coordinates, events arrive. The final background estimate for this 
integration period is the direct convolution of the efficiency map with the rate. 

B{a, S) = J eih, 6) ■ R{a - h)dh (1) 

We refer to this method as ''Direct Integration” since the local event distribution is measured 
and convolved with the detector’s all-sky event rate. The time-independence of e{h, 6) and the 
spatial-independence of R{r) is key because we can use data from the entire integration duration 
in the computation of e and data from the whole sky in the calculation of R, This method has 
been reliably demonstrated to estimate backgrounds with systematic errors of a few parts in 10""^. 
The limiting systematic error is due to real non-uniformities in the cosmic-ray background. 

After computing the background estimate B{a,5), we can take the signal map S{a,6), which 
is just a histogram of arrival directions and compute the excesses by computing S{a.5) — B{a.S) 
in bins of a and 5. 


2.1. A5: Gamma/Hadron Separation Parameter 

We use an event parameter A5 to statistically discriminate air showers induced by gamma rays 
from those induced by hadrons. A5 is defined as 


A5 - 400 ♦ 


MaxPEMU 


The parameter R measures the size of an event and is defined as 


(2) 


Nas Nqr 


( 3 ) 


where Nas/N[1%^' is the fraction of live PMTs in the top layer (or air-shower layer) wTich 
participated in the event and Nq^/N]^ is the fraction of live outriggers to participate in the 
event. The R parameter functions is an estimate of the event’s energy and is described more 
in Section 3. The parameter Ffit is a parameter of the shower-fitting algorithm indicating what 
fraction of the PMTs registered times close to the fitted shower plane. MaxPEMU is the number of 
photo-electrons recorded in the hardest-hit channel from the bottom layer of the experiment. The 
parameter (*(t) is a few percent run-dependent correction to Ffif The distribution is seen to vary 
systematically in the data depending on immodeled factors like changes in the calibration. The 
(^(t) is a correction to take out this variation in Ff^f, The MaxPEMU in the denominator of A5 
is expected to be typically largr^r for hadron- induced air showers because the penetrating particles 
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illuminate the bottom layer of the experiment. This means that Ab is typically larger for gamma- 
ray induced showers with the same number of particles reaching the ground. The numerator of A5 
increases with the size of the event to account for the fact that we expect more light in the bottom 
layer when the event is larger and have to take out the dependence on the overall size of the event. 
The overall scaling factor of 400 gives A5 typical values between 1 and 10. Figure 2 shows the 
distribution of A5 for events in a small circle around the Crab and the separation that A5 provides. 


2.2. Event Weighting 

The Ab parameter provides separation between gamma rays and hadrons primarily because 
of the higher characteristic value for gamma-ray sources. To maximize the statistical significance 
when searching for sources, we assign each event a weight based on its Ab value and the signal-to- 
background expectation for a Crab- like source for events with that Ab. A different set of weights 
is used for each bin. In this approach, more-gamma-like events are counted with a higher weight 
than less-gamma-like events. A hard cut on the gamma/hadron parameter, which is used for the 3 
highest T bins, is simply a step function weight. Weighted sky maps are constructed from data in 
9 T bins between T of 0.2 and 2.0 in steps of 0.2. 

In addition to the Ab weighting for gamma/hadron separation, events are given a weight to 
account for the angular resolution of the instrument. For a given source position hypothesis, an 
additional weight is applied to each event which is a function of the angular distance from the source 
position to the reconstructed event position. We assume a 2-dimensional Gaussian as the form of 
the angular resolution function, where the resolution depends on the event energy parameter, T, 
and ranges from 1.2^ for small T events to 0.35° for large T events. 


2.3. Probability Estimation 

In the absence of weighting, events from signal and background samples are compared and a 
probability for the observed data under the mill hypothesis can be reliably computed using Equation 
17 of Li & Ma (1983). When weighting events rather than simply counting them, we complicate 
the calculation of the expected fluctutions. In the large N limit, this problem has been solved. If 
one records not just the sum of the weights, ™ YLi but also the sum of the squares of the 
weights, N 2 ~ "^'be error in N is computed as 5N — 

In the small N limit however, the SN ~ V^A^2) approximation breaks down, hence the ad- 
vantage of the approach of Li and Ma who derived their probability equation assuming Poisson 
fluctuations. Poisson distributions take discrete values (integers), unlike continuous Gaussian distri- 
butions. Fay & Feuer (1997) point out that the Poisson distributions can be reliably approximated 
as a continuous envelope function described by a single parameter, which for a sum of weights is 
~ Since fluctuations are well approximated as Poisson in , we can rewrite the Li 
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and Ma expression for significance of an observed result as 
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(6) 


and a is the usual ratio of the signal and background exposure. We have studied this approach both 
through examination of data and with Monte Carlo simulations and found it to be reliable even in 
the regime of small statistics, ~ 1. Figure 3 shows the significance distribution for the Milagro 
sky. Since most of the sky has no gamma ray sources, significances are distributed normally. The 
fitted mean between ±2a is -0.013(7 and the width is 0.996(7. This is high-level confirmation that 
the significance calculation is correct. The high-significance tail to this distribution is due to the 
presence of real sources in the sky. Figure 4 shows the final significance map in the region of the 
Crab. The significance at the Crab location is 17.2 standard deviations (tr). This figure includes 
all data over the 8-year operation of Milagro. Only data taken after the outriggers were added are 
used to measure the energy spectrum in this paper. 


3. Energy Estimation 

When a cosmic ray or gamma ray interacts in the atmosphere the amount of energy detected 
at the ground depends on the energy of the primary particle and the depth of the initial interaction. 
Since the Milagro detector is a large- area calorimeter, it is possible to measure the energy reaching 
the ground level with a relatively small error (^20%). However, fluctuations in the longitudinal 
development of air showers - due primarily to fluctuations in the depth of the initial interaction - 
limit the resolution of EAS arrays. Gamma rays of a given energy that penetrate deeply (a few 
radiation lengths) into the atmosphere deliver substantially more energy at the ground level than 
showers of the same energy that that interact at the top of the atmosphere. These fluctuations are 
log-normal (Smith 2008) and dominate the energy resolution for EAS arr^^s such as Milagro. Data 
from September 2005 to March 2008 have been used in determining the energy spectrum because 
the outriggers are needed to provide a dynamic range spanning 1 to 100 TeV. In this dataset, the 
statistical significance of the Crab is 13.5 standard deviations. 
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3.1. The T Parameter 

Figure 5 shows the typical dependence of T , defined in Equation 3, on the primary particle 
energy. We note that a single bin covers a wide range of energies and that these energies overlap 
significantly. Consequently there is no advantage to a finer segmentation than the 9 bins chosen. 

T is well-modeled by the simulation as seen in Figure 6. Shown is the experimentally-measured 
T distribution for background cosmic-rays with the simulation expectation overlaid. The gamma- 
ray enhancing event weights from Section 2.2 have been used as a way to probe the data and 
simulation agreement under the same conditions as eventual gamma-ray spectral measurements. 
Note that the inclusion of the gamma-ray weights significantly restricts the number of simulation 
events surviving to the highest T bins, where the gamma/hadron separation performs the best. 
The weights are, after all, designed to de-emphasize hadronic events. The expected background 
passing rate above of 1.6 cannot be reliably estimated since no simulated background events 
survive the weighting. 


3.2. Spectral Fitting 


Since the energy resolution of Milagro is broad, typically 50%-100%, the energy distribuion 
expected in a bin is dependent on the spectral assumptions. For this reason, we perform all of the 
spectral fits in the T space: For each spectral hypothesis, we simulate the expected T distribution 
and evaluate the goodness of fit based on the statistic. 

We perform spectral fits to a generalized assumption for spectral shape described by a power 
law with an exponential cutoff. 

dN E —E 

-/o(^r"exp(^) (7) 

arj r^Q ^cut 

In this equation, /q, a and Ecut hi parameters for the flux, spectral index and cutoff energy 
respectively. The Eq parameter is not fit but rather is chosen so that the contours of the fit 
variables are de-correlated. This functional form has the benefit that it intrinsically models a pure 
power law hypothesis when Ecut is above a few hundred TeV and we can test a pure power law 
hypothesis and a power law with an exponential cutoff hypothesis with one y^ computation. 

The fit is performed by computing y^ for a given E distribution defined as 


(Jo, a, Ecut ) 


i~jF bins 


a, Ecut ^ Declination) 

“““““ ~5WTsW 


(8) 


Here P and M are the sum of the predicted and measured sum of weighted events per day 
from the Crab and 5P and 5M are the error in P and M respectively. We have only considered 
statistical errors in the estimation of y^. Notice that since the value of P depends on the zenith 
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angle of the source as it transits and the distribution of zenith angles averaged over a transit is 
the same for all sources with the same declination, the predicted daily weight sum depends on the 
declination of the source in addition to the hypothesized spectral parameters. 

The expected weighted excess is computed for discrete values of a between 1.5 and 3.5 in steps 
of 0.025, logio{Ecut) between 0 and 3 in steps of 0.05. Iq was scanned over a range between 0.1 
and 4.5 times the nominal pure power law Crab flux measured by HESS (Aharonian et al. 2006) in 
steps of 0.05. The values are tabulated and the best fit spectrum is computed by minimizing 

Section 4 summarizes the results of this technique applied to the excess from the Crab Nebula 
and to the background cosmic~ray population as a cross-check. 


4. Results 

The success of our technique depends on the simulation to reliably describe the response of the 
instrument. Below 20 TeV, the energy spectrum of the Crab has been well measured by lACTs. 
In the range from 20 TeV to 100 TeV the data are limited and somewhat contradictory. We 
can however test the energy estimation by fitting the spectrum of the hadronic background as a 
cross-check of the method. 


4.1. Systematic Effects 

The spectrum of the hadronic background has been well measured by a series of balloon-based 
spectrometers as well as ground-based air shower detectors. See Particle Data Group et al. (2008) 
for a comprehensive review. While the simulation of hadronic interactions introduces a systematic 
error that is not present in simulated gamma-ray cascades, comparisons with hadronic data are 
nevertheless a useful verification of the simulation. A single day of data is sufficient to fit the 
cosmic-ray background spectrum, so reliable and accurate daily fits to the cosmic-ray spectrum 
serve as a measure of the stability of the energy response of the instrument. 

The hadronic background is composed of numerous species. Protons dominate the flux, ac- 
counting for about 60% of the triggers in Milagro, but helium ions at about 30% and heavier 
element at about 10% also make important contributions. We have simulated the 8 hadronic 
species with the largest contribution: H, He, C, O, Ne, Mg, Si and Eb. The ATIC spectrometer 
has measured both the spectrum and the composition of multi- TeV cosmic rays and found different 
spectra for the different species. We simulate the 8 species listed with their spectra measured by 
ATIC (Panov et al. 2006; Ahn et al. 2006) and fit to an overall offset in the spectral index (Aa) 
and a flux scale factor (S) where Aa™0 indicates that the spectrum was measured to be exactly 
ecpial to the ATIC spectrum, Aa >0 indicates a steeper spectrum and Aa <0 indicates a flatter 
spectrum. Similarly, 5—1 indicates an agreement with the predicted flux, 5 < 1 indicates that 
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Milagro measures a flux that is less than the the predicted flux and 5 > 1 indicates a greater flux 
than predicted. 

Recall that the eventual gamma-ray fits are performed using events weighted by the gamma- 
ray selection weights from Section 2.2. In order that the study of cosmic-ray fits are subjected to 
the same systematic effects that may be present in the gamma-ray analysis, we use the same event 
w^eighting for the cosmic-ray fits, with the consequence that the majority of the cosmic-ray events 
are given small weight and cosmic-ray showers which appear similar to the gamma rays receive the 
most weight. 

Figure 7 shows the fit cosmic ray flux scaling and spectral index as a function of time. Eq was 
chosen to be 10 TeV for these fits. The cosmic ray index varies by less than ±0.1 over the time 
shown. The overall flux scaling changes as the operational conditions of the experiment change. 
Many of these changes are not included in the simulation. Departures from the average are rare 
and suggest a systematic uncertainty in the total flux of sources close to the (Abdo et al. 2007) of 
30% that has been estimated before. 

The instability of the cosmic-ray fit over time is due to real ~ 10% changes in the E distribution 
of the data over time. These changes can be seen easily using the background cosmic-ray data which 
have small statistical errors. Figure 8 summarizes our uncertainty in the E distribution based on 
variations observed in the experimental data. For each of a set of data runs covering the observation 
period, we compute the E distribution of the background data. For each bin of E we quantify the 
width of the distribution of weighted event rates in that bin across the different runs as the 68% 
spread around the median. The fractional width of these distributions for each E are shown as the 
darkest band in Figure 8. Some of the run to run variation is due to an overall scaling difference 
between the runs. If we normalize the E of each run to unit area and re-do the calculation of 
the spread across the runs, we get the darker gray band in Figure 8, It is naturally somewhat 
smaller than the darkest band because variation that can be attributed to overall scaling has been 
taken out. Finally the lightest gray band shows the fluctuations expected due to purely statistical 
effects and we can see that there is a fundamental uncertainty due in the E distribution due 
simply to statistically significant differences in different days of data at the level of about 10%. 
This variation is due to real and unmodeled changes in the detector calibration, configuration and 
operating conditions. 


4.2. Spectrum of the Crab 

Applying the method of Section 3 to the statistical excess from the Crab Nebula, we can 
determine the spectrum of the Crab Nebula. The y* space that is spanned by the three fitting 
parameters from Equation 7 is scanned to find the global minimum. To fit a simple power law or 
to test a number of specific hypotheses motivated by measurements from other experiments, one 
of the parameters can be fixed to the assumed value and the minimum is then found over the 
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corresponding subset of the fit space. 

To begin with, we fit the Crab spectrum under the hypothesis that the spectrum is a pure 
power law. That is to say that the £"cut of Equation 7 is much higher than the Milagro sensitivity. 
The best fit occurs at lo— (6.5 ± 0.4(gtat))^16""^'^(cm^ • s ■ TfeV)“^ and == 3.1 ± O.l(stat) with Eq 
of 10 TeV. For this hypothesis, we obtain a of 24.1 with 7 degrees of freedom. The contours of 
the function are ellipsoidal in the space of the fit parameters indicating very little correlation in 
the fit parameters. Assuming this hypothesis is right, we expect to have only a 0.1% probability to 
observe a by chance. The moderate failure of the two-parameter model to fit the observed data is 
robust even if we artificially inflate the error bars in the data by 10% (added in quadrature) to allow 
for our systematic uncertainty in the rate of events in a given T bin. With the artificially inflated 
error bars, the y^ improves to 21.7 which corresponds to a chance probability of 0.3%. Figure 9 
shows the E distribution for the Crab with our best-fit pure power law hypothesis overlaid. 

An independent analysis of the Milagro data was done (Allen 2007) utilizing a different algo- 
rithm to estimate the gamma ray energy of each event which depended on: the core distance of the 
air shower from the center of the Milagro pond, the reconstructed zenith angle of the air shower 
primary, and the measured number of PMTs in the top layer and outrigger array. Gamma rays 
were distinguished from cosmic rays using the compactness parameter (Atkins et al. 2003) rather 
than A5. The fitted values are consistent with the fits obtained with E and A5 with somewhat 
larger error bars. The agreement indicates that our reported fit is robust with respect to energy 
algorithm and hadron rejection parameter. 

We next consider a hypothesis of a power law, with an exponential cutoff. This is Equation 7 
where Ecut is allowed to vary. With this additional free parameter, the y^ improves to 12,1 with 
6 degrees of freedom. This corresponds to a chance probability of 6%. At the location of the best 
fit, X 10-12 (cm2. s-TeV)-i « = 2.5 ± and TeV. 

For this fit Eq was set to 3 XeV. Figure 11 shows projections of the 1 and 2(j allowed regions in 
the plane of our three fit variables. The somewhat broad allowed range of spectral indices and 
cutoff energies is due to a fundamental ambiguity in the Milagro data that a soft spectrum is hard 
to distinguish from a harder spectrum with an exponential cutoff. Fixing the low-energy spectral 
index to the values between 2.4 to 2.6, as measured by other experiments, gives a la allowed range 
for the cutoff energy of between 20 and 40 TeVb 

Neither of the two spectral assumptions is preferred strongly by fitting the Milagro data. The 
pure power lavr fit is a mmginally poor fit. The addition of an exponential cutoff improves the 
fit. The measured fluxes are shown on Figure 10 for the two hypotheses. Regardless of which fit 
is chosen, the general conclusion is clear: The high-energy spectrum, above about 5 - 10 TeV, is 
steeper than measured by lACTs at lower energies. In the pure power-law hypothesis, this manifests 
itself as a measured spectrum of a ~ 3.1 ± 0 . 1 , steeper than has been measured by, for instance, 
HESS of 2.4 to 2.6. The fit that allows for an exponential cutoff reproduces the low-energy spectral 
index measured by LACTs and this steepening at high energy is seen as an exponential cutoff at ~ 
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30 TeV. 

Finally, it is interesting to note that above 30 TeV, the HESS and HEGRA data are mildly 
inconsistent. The HEGRA measurement continues to higher energy than the HESS data. It has 
been suggested (Bednarek & Idee 2011) that this discrepancy is related to the time variability 
observed by the Crab since HEGRA was observing earlier than HESS. The Milagro data, which 
represents the time-average over 3 years of data, indicates a spectrum between the data of HESS 
and HEGRA. 


5. Conclusions 

The Crab Nebula is the brightest northern hemisphere TeV source and has been extensively 
measured by lACTs above 1 TeV. The Milagro measurement of the energy spectrum of the Crab 
has been presented. A background rejection parameter (A5) has been described and shown to 
distinguish between gamma ray and hadronic primaries in the detector. We have presented the 
weighting and background estimation and background subtraction techniques used to extract the 
Crab signal, giving a 17a over the lifetime of the experiment. 

The size of an air shower at ground represented by the fraction of PMTs in the Milagro 
experiment that detect a signal (the T parameter), is a suitable variable for measuring the spectra 
of primary TeV gamma and cosmic rays. The relatively simple form of T is justified because the 
dominant effect contributing to the energy resolution of Milagro is fluctuations in the depth of 
first interaction of the primary particles and not in the measurement of the energy reaching the 
ground. The parameter is well modeled in the simulation as observed by studying the cosmic ray 
background. 

The energy spectrum of gamma rays from the Crab between 1 and 100 TeV has been measured 
by fitting the observed T distribution of the Crab with expectations from simulation. A steepening 
of the spectrum above about 5-10 TeV with respect to measurements by lACTs at lower energies 
has been measured. 

The experiment observes the entire overhead sky, the data and analysis technique presented 
here for the Crab observations can be used to measure the flux and spectral properties of the other 
sources in the Milagro catalog. The agreement seen on the Crab as a calibration source justifies 
confidence in measurements of other sources. 
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Fig. 1. — Example of the background subtraction technique in a single declination bai 
display purposes, this calculation is performed with a four-hour integration instead of the s 
two-hour integration. R{t) is the all-sky event rate. e(h,S) is the local-coordinate distrib 
event mrivals and is convolved with R(t) to arrive at B{a,S), the background estimate, B 
subtracted from the binned event arrival directions S(a. S) to arrive at the actual excess esi 
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Fig. 2. — Shown is the distribution of A5 compared to simulation, both for the background cosmic 
rays and for the background-subtracted Crab excess. The gamma ray signal is shown in a small 0.7 
degree circle around the true Crab location. The A5 parameter is used to distinguish gamma-ray 
events from hadronic events. A higher value of A5 indicates a higher probability than an event 
originated from a gamma ray. 
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Fig. 4. — Shown is the statistical significance in the region around the Crab Nebula, indicated by 
the white dot. The gamma-ray-enhancing weights as well as the angular smoothing from the text 
have been incorporated. Data over the entire 8-year lifetime of the experiment have been used and 
all T bins have been combined. At each point in the map, the statistical significance is calculated. 
The smoothing causes the points to be very correlated. The significance at the location of the Crab 
is 17.2tr. 
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Fig. 5. — Typical dependence of J' with eiF 
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ys A source with a spectrum of 2.6 is assumed. Shown 
orgies for events in the indicated F range. 
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Fig. 6. — The distribution of the background cosmic rays in a 2.5 degree bin at the same 
declination as the Crab, but offset by a few degrees in right ascension. Note that the weights from 
Section 2.2 have been applied in this comparison as a way to probe potential systematic biases 
introduced by the weighting procedure. The statistical errors in the simulation expectation for 
high F are quite large due to the weighting which leaves very few cosmic ray events with high 
weight in highest T bins. 
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Fig. 7, — The left figure shows the fitted cosmic ray scale (relative to the simulation) as a function 
of time and the right figure shows the fitted cosmic ray index, as an overall steeping or flattening 
of the simulated spectrum. Note that T greater than 1.4 was not used in these fits. 
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Fig. 8. — This figure quantifies uncertainty in the T distribution. The 68% central values for back- 
ground weighted event rates going into JF bins across different days of data acquisition are shown. 
Both the absolute variations are shown along with variations after normalizing the underlying T 
distributions to unit area. This indicates a fundamental 10% uncertainty in the shape of the F 
distribution. 
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Fig. 9. — The distribution of the T parameter measured in data and expected from simulation for 
several hypotheses. The first two are the best fits to the Milagro data both with a pure power-law 
and power-law with an exponential cutoff. The second two show the expectations from the best 
HESS solutions, both for a pure power law (with a differential photon spectral index of -2.63) 
and including an exponential cutoff (with an index of -2.39 and a cutoff at 14.3 TeV) as reported 
in Aharonian et al. (2006). Note the points have been offset horizontally by a small amount for 
display. 
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Fig. 10. — The panels display the spectrum of the Crab as measured by R4ilagro. The pure power 
law hypothesis is shown on the left and the power law with an exponential cutoff is shown on the 
right. The one a regions ai'e shown with shading and the points measured by HESS, HEGRA and 
Whipple are overlaid. 
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Fig. 11. — Shown are the 1 and 'la allowed regions for the full power- law hypothesis including an 
exponential cutoff. The position of minimum is indicated. The 3-d regions have been projected 
into planes of the three fit variables. 







